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Abstract 

The existence and stability of solitons in Bose-Einstein condensates with attractive inter-atomic 
interactions, described by the Gross-Pitaevskii equation with a three-dimensional (3D) periodic 
potential, are investigated in a systematic form. We find a one-parameter family of stable 3D 
solitons in a certain interval of values of their norm, provided that the strength of the potential 
exceeds a threshold value. The minimum number of Li atoms in the stable solitons is 60, and 
the energy of the soliton at the stability threshold is ~ 6 recoil energies in the lattice. The 
respective energy-vs.-norm diagram features two cuspidal points, resulting in a typical swallowtail 
pattern, which is a generic feature of 3D solitons supported by low- (2D) or fully-dimensional lattice 
potentials. 

PACS numbers: 03.75.Lm,03.75.Kk,05.45.Yv 
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Creation of multidimensional solitons built of nonlinear light or matter waves is a great 
challenge to the experiment. The current situation in this field is summarized in a recent 
review The only real example of a quasi- two-dimensional (quasi-2D) spatiotemporal 
soliton in optics was one created in a (quadratically nonlinear) crystal |2|; it could not 
be made fully three-dimensional, as one transverse direction was reserved to implement the 
technique of tilted wave fronts, by means of which sufficiently strong artificial dispersion was 
induced in the medium. 

While the experimental situation in optics remains difficult [|, new possibilities for the 
creation of multidimensional solitons are offered by Bose-Einstein condensates (BECs) with 
attractive interactions between atoms. In an effectively ID condensate of 7 Li, single solitons 
and soliton clusters were successfully created |3jj. In the 2D and 3D situations, the self- 
focusing cubic nonlinearity induced by the attractive interactions leads, respectively, to weak 
and strong collapse, as predicted by the corresponding Gross-Pitaevskii equation (GPE) 
[alias the nonlinear Schrodinger (NLS) equation] 4], and demonstrated experimentally in 
the BEC 0. However, a spatially periodic potential in the form of an optical lattice (OL), 
created as an interference pattern by coherent laser beams illuminating the condensate, can 
stabilize the multidimensional solitons in the self-attractive BEC. For the 2D case, this was 
& ,t predicted in Ref, RQ; _, it was also debated Q that 2D solitons may 
be stable in the presence of a quasi-lD periodic potential (one that does not depend on 
;he second spatial coordinate). 3D solitons supported by the 3D OL were reported too 
f|. The form of the soliton was predicted by the variational approximation (VA), which 
was used, as an initial guess, to generate several examples of stable 3D solitons in direct 
simulations, including "single-cell" and "multi-cell" ones, that are essentially confined to 
one or several cells of the OL, respectively. Further, in Refs. 8| and it was independently 
demonstrated that stable 3D solitons can also be supported by the low-dimensional, i.e., 
quasi-2D, lattice potential (however, the quasi-lD potential cannot stabilize 3D solitons jsj], 
unless it is combined with periodic alternation of the nonlinearity sign, provided by the 
Feshbach resonance in ac magnetic field ^^). In particular, the VA predicts that the 3D 
solitons in the quasi-2D lattice exist at all values of the norm (number of atoms) N and 
OL strength p, but are stable only for p exceeding a certain critical value p CI , in an interval 



of the width AA^ 3 ^ ~ (p — p CT ) 1 ^ 2 in a vicinity of a finite value N Cf of the norm |8[ [the 
stability was predicted on the basis of the Vakhitov-Kolokolov (VK) |Tj| criterion]. 
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Except for a few examples reported in Ref. |6j, no systematic investigation of the ex- 
istence, stability and robustness of 3D BEC solitons in the 3D lattice potential has been 
performed yet. The objective of the present work is to report basic results for this problem. 

The GPE provides for the description of the BEC dynamics in terms of the mean-field 



single-atom wave function ip(x,y, z,t) [12[. The normalized form of this equation for a 
self-attractive condensate trapped in the 3D potential —V(x,y,z) is well known 



dt 2 \dx 2 dy 2 dz 2 
Generally, V(x, y, z) contains terms accounting for the confining parabolic trap (magnetic 

and/or optical) and the periodic potential of the OL. Being interested in the localized 
solutions, occupying a few cells of the lattice, we disregard the parabolic potential and 
set V(x, y,z) = p[cos(4x) + cos(%) + cos(4z)], where the OL period is normalized to be 
7r/2, and the OL strength p is defined to be positive. Besides the above-mentioned norm, 
N — J J J \ip(x,y,z)\ 2 dxdydz, Eq. (JTJ conserves the energy E (see [l^)- Stationary soli- 
ton solutions have the form tp(x,y, z,t) = w(x, y, z) exp(— ifd), with a real function w and 
chemical potential \l. We looked for the function w(x, y, z) by means of the known method 
of the propagation in imaginary time ^j]. It was implemented, using the Crank-Nicholson 
scheme, with the nonlinear finite-difference equations solved by means of the Picard iteration 
method, and the resulting linear system handled with the help of the Gauss-Seidel iterative 
procedure. To achieve good convergence, we typically needed six Picard iterations and six 
Gauss-Seidel iterations. We used equal transverse grid stepsizes, Ax = Ay = Az = h, and a 
mesh of 361 x 361 x 361 points was usually employed. The convergence to a stationary state 
occurred after 4 x 10 3 — 5 x 10 4 steps of the evolution in imaginary time, typical transverse- 
grid and time stepsizes being h = 0.02 and At = 0.0003, respectively, for narrow solitons, 
whereas for broad ones it was enough to take h = 0.07 and At = 0.004. 

One can derive a relationship between the norm N, chemical potential /i, real wave 
function w and energy E of the stationary solution: E = fiN + | J f J w A (x,y, z)dxdydz. 
It can be used to determine the chemical potential fi, once the field profile w is known. 
This exact relation was also used for verification of accuracy of the numerically found sta- 
tionary solutions. Notice that, for stationary solitons of the NLS equation in the free 3D 
space [with V = in Eq. JJJ], the following relations between /z, N, and E are known: 
fi(N) = -CN- 2 , E(N) = CN~\ C w 44.5 Q,Q (a corollary of this is dfi/dN > 0, which 



immediately shows that these free-space solitons are always unstable, as the VK stability 
criterion requires exactly the opposite, dfi/dN < 0). 

In Figs. 1(a) and 1(b) we plot the dependences /i = fi(N) and E = E(N) for the 
numerically found family of 3D solitons in the present model. It is seen that, in the presence 
of the 3D OL, the localized states exist only for [i smaller than some maximum value, /i max (p)j 
in fact, it corresponds to the edge of the bandgap in the spectrum of the linearized equation 
(JTJ). At values of fi that do not belong to the bandgap, no soliton is possible [l^. We note 
that /i max (p) decreases with the increase of the lattice strength p, see Fig 1(a) [for the 3D 
NLS equation in free space (p = 0), one has /i m ax(0) = 0]. Remarkably, for sufficiently large 
values of the lattice strength p, the E(N) curves in Fig. 1(b) feature two cusps, instead of a 
single one, as in most other 2D and 3D Hamiltonian models. Examples of the usefulness of 
the energy-vs.-norm diagrams (or Hamiltonian-vs.-power diagrams, in the context of spatial 
optical solitons) in the analysis of the existence and stability of solitons can be found in Ref. 
|l6j | . To check the soliton's stability by direct propagation simulations we used the standard 
Crank-Nicholson discretization scheme with 121 x 121 x 121 points in the coordinates x, y, z, 
and spatial and time stepsizes 0.019 < h < 0.073 and 0.00045 < At < 0.004 (depending 
on the soliton's size). The validity of the VK stability criterium was checked by performing 
direct propagation simulations; the step in the soliton's parameter was A/i = 0.08. Thus, 
we have found stable 3D solitons, in a finite interval of the values of N, if the OL strength 
p exceeds a threshold value p cr : as seen from Fig. 1, p cr itself is located between p — 1.5 
and p = 2. Further, Fig. 2 displays an integrated characteristic of the soliton family - 
the dependence E = E(fi, N) - together with the corresponding projections onto the three 
planes (E,fi), (E,N), and (fi,N), for two different values of the lattice strength, p = 
and p = 3. In the latter case [in Fig. 2(b)], we notice the non-monotonic behavior of the 
three projected curves for p = 3 (above the stability threshold) and the "swallowtail" loop in 
the energy-vs.-norm diagram. Although this pattern is one of generic possibilities known in 
the catastrophe theory, it rarely occurs in physical models (applications of the catastrophe 
theory to the soliton-stability problem were reviewed in Ref. (llj). 

It is noteworthy that the qualitative features of the stability picture for the 3D solitons are 
essentially the same as found earlier in the approximate analytical form by means of the VA 
11 , and in the numerical form as well 0,0], for 3D solitons supported by the low-dimensional 
(quasi-2D) lattice: a narrow stability interval AN appears in a vicinity of a finite critical 
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value N CI , when the lattice strength p exceeds a finite minimum value p cv . Remarkably, 
the "swallowtail" loop was also identified as a characteristic feature of the family of stable 
3D solitons supported by 2D harmonic lattices jjj. It is relevant to compare the stability 
picture for the 3D solitons supported by a fully- dimensional OL with that for 2D solitons 
supported by an OL (that may be either a low- dimensional quasi-lD lattice, or the full 2D 
one) 0, Is, la]: in that case stable solitons appear at arbitrarily small values of the lattice 
strength, and their stability interval extends, in terms of N, up to a maximum ("cutoff") 
value at which the 2D solitons exist (the latter is actually the norm of the Townes soliton 

n 

i.e., a soliton of the free-space 2D NLS equation [4]). Thus, we conclude that all the above 
features, viz., the existence of stability threshold p CT , the (small) stability interval AiV sta b for 
p > p CT , and the typical swallowtail pattern in the dependence E = E(N), are generic to 3D 
solitons supported by lattice potentials, and distinguish them from 2D solitons. 

Shapes of both unstable and stable solitons are shown in Fig. 3, through their isosurface 
plots for a typical value of the lattice strength parameter, p = 3. An unstable low-amplitude 
3D soliton (with A = w(0, 0, 0) = 1.8), found for the value of the chemical potential fi close 
to the bandgap edge, is displayed in Fig. 3(a). This soliton is broad, occupying many lattice 
cells. Typical stable solitons, with medium (A = 2.2) and high (A = 3) amplitude, are 
shown in Figs. 3(b) and Fig. 3(c), respectively. Notice that the unstable soliton in Fig. 
3(a) and its stable counterpart in Fig. 3(c) have equal values of the norm, N = 2.4. The 
intermediate stable soliton in Fig. 3(b) has a smaller norm (N = 2.04), which is very close 
to the limit value corresponding to the first cuspidal point in the dependence E = E(N), 
see the curve pertaining to p = 3 in Fig. 1(b). Figure 4 additionally shows integrated views 
along the z-axis of the isosurface plots displayed in Fig. 3. This figure illustrates the fact 
that low-amplitude solitons spread to many lattice cells, whereas high-amplitude solitons 
occupy only a few cells. 

An important issue for these 3D solitons is occurrence of the collapse (recall that the 
3D collapse in the NLS/GPE is strong, on the contrary to the weak 2D collapse We 
expect that solitons on the stable branch are able to withstand small perturbations without 
collapse, whereas linearly unstable solitons would either collapse, or reshape themselves into 
time-periodic breathers, or decay into "radiation" (quasi-linear Bloch waves), depending 
on the type and strength of the perturbation. In order to verify these expectations, we 
simulated the evolution of the solitons under small perturbations, taking the initial condition 
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as ij)(t = 0) = w(x, y, z)(l + ep), where e is a small amplitude of the perturbation, and p 
is taken either as a random number from the interval [—0.5, 0.5] (stochastic perturbation) 
or as p = 1 (uniform perturbation). We have checked that the solitons belonging to the 
VK-stable segments of the curve p = p{N) are indeed stable against small perturbations. 
To illustrate this result, Fig. 5 shows an example of a stable soliton which persists after the 
application of the stochastic perturbation with e = 0.1: in the course of the evolution, the 
soliton's amplitude slightly oscillates, with no trend to collapse or breakup. This and many 
other simulations clearly show that the linearly stable solitons are also nonlinearly stable 
objects. For linearly unstable solitons, the simulations reveal the following scenarios of the 
instability development, (i) Low- amplitude solitons decay into linear waves under stochastic 
perturbations, or under uniform ones that reduce the soliton's norm, i.e., perturbations with 
p = 1 and e < 0; this case is illustrated by Fig. 6 for a typical situation, (ii) The same 
unstable soliton, but under uniform perturbations with p = 1 and e > (that make its 
norm larger), reshapes itself into a time-periodic breather, (iii) An unstable high-amplitude 
soliton collapses if it is perturbed by the uniform perturbation that increases its norm. 

In conclusion, we have found one-parameter soliton families of the 3D nonlinear 
Schrodinger/Gross-Pitaevskii equation with self-focusing nonlinearity and 3D lattice po- 
tential, and explored their stability. Comparing the results with recently published find- 
ings for stable 3D solitons in the quasi-2D lattice potential, and with (very different) re- 
sults for 2D solitons supported by the quasi- ID or full 2D lattice, we were able to iden- 
tify generic features that distinguish the stable 3D solitons: a finite stability threshold 
p cr in terms of the lattice strength p, a finite stability interval AiV in terms of the soli- 
ton's norm N, and the swallowtail shape of the energy-vs.-norm dependence. The 3D 
solitons investigated here can be created in BEC - most plausibly, using 7 Li loaded into 
a 3D optical lattice, with the spatial period A ~ 0.5 pm. Then, undoing normalizations 
that cast the GPE into the rescaled form it is easy to find that the actual number 
of atoms iVphys is related to the solution norm by iVphys = (27r 2 |a|) 1 AiV, where a is the 
scattering length of the atomic collisions. For 7 Li, a = —1.45 nm, which yields the mini- 
mum number of atoms necessary for the formation of stable solitons, that corresponds to 
N CT ~ 3.5 in Fig. 1: ps 60 (note that it corresponds to the density ~ 5 x 10 14 

cm -3 , which is quite a typical value for BEC). The strength e of the OL, chemical poten- 
tial, and energy in physical units are related to their normalized strength counterparts as 
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e = (l/8)E rec p, /i phys = (l/4)£ rec //, E phys = [H 2 / (8\a\Am)]E = {A/[(47r) 2 |a|]}£ rcc £, where 
E TCC = (2 /to) (nh/ A) 2 « 7.5 x 10~ 29 J is the recoil energy in the lattice. Thus, near the 
soliton-stability threshold, we get e ~ 0.25i? rec , /i p h ys ~ —E rcc , and -Ephys ~ — S-E'rec- 
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FIG. 1: The chemical potential \i (a) and the energy E (b) of the 3D solitons versus their norm 
TV for different values of the lattice strength p. Full and dotted lines depict stable and unstable 
solitons, respectively. 

FIG. 2: (Color online) The soliton family in terms of the dependence E = E(/j,,N), and its 
projections on the planes (E, fi), (E, N), and (//, N) for p = (a) and p = 3 (b). 

FIG. 4: (Color online) Integrated through views along the z-axis of the solitons shown in Fig. 3. 



FIG. 5: Isosurface plots showing self-cleaning of a stable soliton corresponding to p = 3 and 
N = 2.4, initially perturbed by white noise, (a) Input at t = 0; (b) output at t = 50. 



FIG. 6: (Color online) Integrated through views along the z-axis of an unstable soliton that decays 
due to a uniform norm-reducing perturbation with e = 0.01, for p = 3. (a) t = 0, A = 1.78; (b) 
t = 50, A = 0.6505,(c) t = 70, A = 0.2075. 



FIG. 3: Isosurface plots of unstable (a) and stable (b) and (c) 3D solitons. Here p = 3; (a) N = 2.4, 
A = 1.8, (b) N = 2.04, A = 2.2, and (c) N = 2.4, A = 3. 
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